The offshore environment, particularly in the deep water case, offers a combination of unique circumstances that make the application of exploration methods based upon remote measurements of electrical resistivity a useful and cost effective exploration tool. These circumstances are:                1. The high cost of offshore exploration wells.        2. The very low level of electrical noise present at the ocean bottom.        3. The large resistivity contrast in many locations between oil filled reservoir rock and the surrounding brine filled sandstone and shale.        4. The excellent electrical coupling provided by seawater.        5. The operational feasibility of collecting resistivity measurements over large area.        
Resistivity is a measure of the amount of electrical resistance offered by a unit volume of material; it is measured in units of ohm-m. Marine resistivity measurements are made by towing a transmitting bipole antenna (100 to 500 m long) with grounded ends approximately 50 m above the ocean bottom and recording the electric field at grounded detecting bipole antennas (typically about 5 m long) mounted on sea bottom recording devices. When a human operated EM source is used, the survey is called a controlled source electromagnetic (“CSEM”) survey. The recording devices measure the electric fields along the two horizontal directions and sometimes along the vertical direction as well. The magnetic fields can be measured as well; in which case, coils are used pointing along the horizontal and vertical directions. The detecting devices store the recorded electromagnetic fields on internal disk drives. When the survey is complete, the recording devices detach themselves from their anchors in response to a signal from the survey vessel and float to the ocean surface. Eventually they are collected by the survey vessel and the recorded data is loaded into onboard computers. The recorded data is correlated with the transmitted signal, transmitter location and heading and other information based upon very accurate independent time recordings (clock readings) on the detector and onboard the ship. Clock readings are used because there is only very minimal communication between the survey vessel and the sea bottom detectors.
Data are usually analyzed in the frequency domain because only certain frequencies are useful due to the skin depth effect which predicts the decay of electromagnetic radiation in conductive materials with propagation distance. The skin depth is the distance over which the radiation carried by a plane wave (a waveform in which the surfaces of constant phase and amplitude are infinite planes perpendicular to the direction of wave propagation) in a uniform material will decay by one factor of e≈2.718 (the base of the natural logarithm system). The skin depth, δ, (expressed in meters) at frequency f (in Hertz), for a material with resistivity, ρ, (in ohm-m) is approximately,
  δ  =                    (                  1          ⁢                      /                    ⁢          2          ⁢                                          ⁢          π                )            ⁢                        (                                    10              7                        ⁢                          p              /              f                                )                      ≈          503      ⁢                        (                      p            /            f                    )                    
In plain words this equation states that the skin depth is approximately 503 times the square-root of the ratio of the resistivity to the frequency. Thus the skin depth will be 503 m in a material of 1 ohm-m resistivity for radiation at 1 Hertz. A target located 1 km below the sea bottom in a 1 ohm-m material will only be illuminated by radiation in the low frequency range ( 1/10 to 1 Hertz). Above 1 Hz, the radiation will experience too much decay in propagating down to the target and back. Below approximately 1/10 Hz the radiation will average too much of the material to effectively sense typical target geometries of interest which tend to be 1 to 10 km by 1 to 10 km in the horizontal dimensions and 10 to 100 m thick. Because the data are analyzed in the frequency domain, the transmitting antenna broadcasts time domain signals which are selected to be rich in only the frequencies of interest. A square wave is often used because it is easy to produce electronically. The period (in seconds) is usually selected to be T=1/f where f is the lowest frequency of interest. The higher harmonic components of the square wave (3f; 5f, etc) help provide resolution of the resistivity variations between the sea bottom and the target.
Time domain solutions to CSEM exist but are computationally even more costly than frequency domain methods. Although they offer clear advantages in handling the so-called air wave effect in land or shallow water surveys (the air wave is the direct transmission from the broadcasting antenna to the detecting antenna through the air), time domain analysis in the deep water case loses the advantages of the frequency domain method because the time domain response is spread across frequencies of little interest to hydrocarbon exploration for targets at typical reservoir depths.
Time domain electromagnetic data (i.e., as measured in a survey), and frequency domain electromagnetic data are related by the Fourier transformation:
      D    ⁡          (      f      )        =            (              1        ⁢                  /                ⁢        2        ⁢                                  ⁢        π            )        ⁢          ∫                        ⅇ                                    -              2                        ⁢            π            ⁢                                                  ⁢            ift                          ⁢                  D          ⁡                      (            t            )                          ⁢                  ⅆ          t                    Because i=√{square root over ((−1))}, it can be seen that electric and magnetic fields in the frequency domain are complex quantities (numbers that have real and imaginary parts). These quantities display characteristic variations with transmitter-to-detector offset, frequency, and transmitter direction which, under appropriate conditions, can be used to constrain subsurface variations in resistivity. The analysis and interpretation of electromagnetic data are carried out by a process of modeling likely geological situations and comparing these simulated data with actual measured data. Very often, the geologic models are adjusted and simulated again and again until an acceptable match between measured and computed data is achieved. When this process is carried out within a computer driven optimization scheme, it is referred to as an inversion.
The interpretation process, either a by hand comparison of models and measured data, or an actual inversion, benefits from the unusually high level of precision with which actual field recordings can be made. Experience now confirms that field recording can be accurate across a dynamic range of approximately seven orders of magnitude. This extraordinary dynamic range makes the CSEM technique sensitive to even the relatively thin targets of interest in hydrocarbon exploration. Unfortunately, the extreme accuracy of the data also makes simulation of the electromagnetic fields a very demanding computational task. Geologic bodies of actual interest are very rarely well approximated by a one dimensional model (layered models) or by a two dimensional model in which the geology is assumed to be unchanging in the perpendicular dimension. Very often the full complexity of a three dimensional (3-D) model must be confronted.
Models of actual interest can be treated by decomposing the 3-D resistivity volume into rectangular cells which are then used to numerically solve the system of equations which govern the electromagnetic fields. Two general classes of such models exist: isotropic models in which the cell conductivities (the reciprocal of resistivity) are independent of current direction, and anisotropic models in which the conductivities depend upon current direction. In general conductivity is a rank 2 tensor with components representing the current along the three Cartesian directions resulting from an applied electric field along each of three directions. The equations relating the electric and magnetic fields to imposed electric and magnetic current sources are known as Maxwell's equations and for typical CSEM data processing application they are expressed in the frequency domain. This approach is called the frequency domain finite difference (FDFD) method.
Two alternatives to the finite difference method (“FDM”) of numerical solution exist for general 3-D models—the finite element method, hereinafter “FEM”, and the integral equation method, hereinafter “IE.” The FEM in which the resistivity grid is composed of more general geometric shapes than rectangles is more complex computationally than the FDM. The FEM is definitely preferred for engineering applications in which great detail is required. Aside from the water bottom topography, the attention to detail offered by the FEM appears not to be required in the application of CSEM to hydrocarbon exploration. Nearly all underlying physical and mathematical details of the FDM also apply to the FEM. The IE method recasts the system of differential equations implied by Maxwell's equations into an associated integral equation by making use of the properties of the Green's function for the electric and/or magnetic field in a uniform or layered model. A uniform or layered material is used for the reference Green's function because highly accurate and rapidly computed solutions are available for these models. The resulting integral equations naturally give rise to computational schemes that work very well for compact 3-D objects imbedded in uniform or layered models (such as a ship or aircraft in the deep ocean or high in the atmosphere). Unfortunately, IE methods appear to be poorly suited for modeling applications when the models exhibit significant variations nearly everywhere, as is expected in the hydrocarbon exploration case. An additional disadvantage of the IE method is that direct discretization of the integral equations gives rise to systems of equations which are highly non-sparse (the electric field unknowns obey equations for which all terms contribute).
Unfortunately, the equations of the FDFD method need to be solved for each frequency of interest and for each of the many thousands of transmitter positions occupied by the broadcasting antenna used in an actual survey. In three dimensions this is an enormous undertaking because of the volume of the earth affected by the broadcast signal. Typical models of interest contain one to ten million resistivity cells when account is taken of the necessary sampling requirements needed to obtain accurate numerical simulations and to avoid contamination of the solution by undesirable boundary effects caused by imposing a resistivity grid of finite size. The computational magnitude presented by a full inversion of an actual survey is so large as to be beyond the capabilities of even the very largest parallel computers, composed of tens of thousands of processors, operating for many months unless computational efficiencies are developed and applied.
One often suggested computational efficiency approach involves recognition of the fact that the transmitter locations are closely spaced (50 to 100 m) so that considerable computational efficiency might be found by simply making small corrections to an already computed solution for a nearby transmitter location. While the possibility exists that some method along these lines might be made to work, the present inventors have experimented several times with this procedure and found that correcting solutions from nearby transmitter locations is just as expensive as obtaining independent solutions for all transmitter locations when seeking solutions that reflect the high accuracy observed in the actual measured data.
Another computational efficiency approach [1] involves generating a complete factorization or decomposition of the Maxwell operator. In a sense this approach amounts to obtaining a solution for all source positions in a form which can be rapidly evaluated for any specific source location. To explain this, some further detail concerning numerical aspects of the FDFD method are required. Maxwell's equations in the frequency domain relate the time-harmonic electric and magnetic fields to the current densities of the impressed electric and magnetic transmitters. As long as one is not concerned with the electrostatic or magnetostatic (the DC) portion of these fields, one may mathematically eliminate the magnetic fields from the equations (this elimination gets into trouble at very low frequency where the electrostatic and magnetostatic contributions are no longer coupled; this is called the null-coupled problem). The resulting equation relates the electric field to the impressed electric transmitters (we assume that our earth model contains no variations in the magnetic permeability which is a good assumption unless ferric materials are nearby) by a system of equations which become a system of algebraic equations for the electric field on the resistivity grid. This system of equations takes the form:Ke=jwhere K is conceptually a square complex matrix of dimension 3N×3N where N is the number of cells in the 3-D resistivity grid (about one to ten million) and e and j are 3N-long column vectors representing the electric field unknowns at the grid nodes and the impressed transmitter electric current. The matrix K is symmetric because of the Reciprocity principle. Fortunately the matrix K is extremely sparse (in fact all but 13×3N of its 3N×3N elements are zero). The above equation is a system of 3N equations for the 3N electric field values on the grid. A fundamental theorem of linear algebra [1] requires that any m×n complex matrix is unitarily equivalent to a diagonal matrix,K=UDV*. 
In this equation, U and V are unitary matrices and D is diagonal. The asterisk denotes the adjoint or complex conjugate-transpose of the matrix V. The transpose of a matrix is a copy of the original matrix in which row and columns have been exchanged,VijT=Vji The complex conjugate of a complex number z=x+iy is the complex number x−iy. The superscript T denotes the transpose and i and j refer to specific matrix elements (i and j assume integer values from 1 to 3N). A diagonal matrix is zero except along the main diagonal. For the special case where K is square, U=V. This is called the singular-value decomposition of K. The non-zero elements of the diagonal matrix D are the singular values of K. Since a unitary matrix is, by definition, one in which the matrix inverse is equal to its adjoint (the matrix inverse is the matrix which produces the unit matrix when multiplied by the original matrix (from either side for square matrices)),U−1=U*
A singular-value decomposition or UDU* factorization of K produces the desired electric field solution via the equation,e=UD−1U*j The singular-value decomposition in effect provides a unitary transformation to a new coordinate system in which the matrix K is diagonal. Once one knows what this coordinate system is, the process of determining the electric field solution for any transmitter is relatively inexpensive. Unfortunately, the computational cost of such decompositions is of order (3N)3 making this approach too costly for 3-D models, but possibly affordable for 2-D models for which N will be approximately 100 times less than in 3-D (and therefore one million times less expensive for a decomposition because the cost is cubic in N). A variant of this approach which involves finding all of the singular values of K is to decompose the matrix K only up to given accuracy by means of an iterative procedure. The singular values in D can be sorted by magnitude from the biggest to the smallest. In fact, many are very close to zero. There exist iterative procedures for producing the decompositions with a cost of (3N)2 operations per recovered singular value. Experiments with using nearby transmitter solutions as the starting point for iterative solutions demonstrate that the effective spectral radius (the range of significant singular values in D) required for accurate electric field solutions may be extremely broad so that this approximate approach may not be cost effective. Again in two dimensions this approach might work. This approach might be effective for the well logging problem because sources are always activated in a restricted portion of the model (the borehole). A further difficulty in the complete decomposition approach is that the most widely used mathematical algorithms do not lend themselves to a parallel implementation.
The present invention provides an alternative approach to efficient processing of 3D CSEM data.